Application of the Space-Time Conservation 
Element and Solution Element Method 
to Shock- Tube Problem 


Xiao- Yen Wang, Chuen-Yen Chow 
Department of Aerospace Engineering Sciences 
University of Colorado at Boulder 
Boulder, Colorado 80309-0429 
and 

Sin-Chung Chang 

National Aeronautics and Space Adminitration 
Lewis Research Center 
Cleveland, Ohio 44135 

Abstract 

An Euler solver based on the method of space-time conservation element and 
solution element is used in this paper to simulate shock-tube flows involving shock 
waves, contact discontinuities, expansion waves and their intersections. Seven test 
problems are considered to examine the capability of this method. The numerical 
results, when compared with exact solutions and/or numerical solutions by other 
methods, indicate that the present method can accurately resolve strong shock and 
contact discontinuities without using any ad hoc techniques which are used only at 
the neighborhood of a discontinuity. 
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Introduction 

The method of space-time conservation element and solution element ( to be ab- 
breviated later as the CE/SE method) is a numerical method developed recently 
by Chang[l] for solving conservation laws governing the motion of fluids. This 
new method differs substantially in both concept and methodology from the well- 
established traditional methods such as the finite difference, finite volume, finite el- 
ement and spectral methods. It is conceived and designed to overcome some of the 
major shortcomings of the traditional methods. 

The development of this new method is guided by several basic requirements. They 
are: (i) to enforce both local and global flux conservation in space and time with flux 

l 

evaluation at an interface of two conservation elements being an integral part of the 
solution procedure and requiring no interpolation or extrapolation; (ii) to use local 
discrete variables such that the set of variables in any one of the numerical equations 
to be solved is associated with a set of immediately neighboring solution elements; 
(iii) space and time are unified and treated on the same footing; (iv) mesh values of 
dependent variables and their derivatives are considered as independent variables to 
be solved simultaneously; (v) to minimize numerical diffusion, a numerical analogue 
should be constructed, as much as possible, to be compatible with the space-time 
invariant properties of the corresponding physical equations; and (vi) to exclude the 
use of the characteristics-based techniques, and to avoid the use of ad hoc techniques 
as much as possible. It was shown in [1] that the above requirments can be met with 
a simple unified numerical framework. 

One of the numerical schemes developed in [1] is an Euler solver for solving the one- 
dimensional unsteady Euler equations. It was used to obtain accurate flow solutions 
involving an infinitely long shock tube. 

In this work, the above Euler solver is extended and applied to more complex 
flow problems involving shock tubes of finite or infinite length. Shock-tube flows 
represent a class of particularly interesting test cases that have exact solutions in 
the presence of shock waves, contact discontinuities, expansion fans and even their 
reflection and intersection when the tube length is finite. The ability to generate 
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numerical algorithms which allow a high resolution of discontinuities, such as shock 
waves and contact discontinuities, without numerical oscillations is an important 
requirement of a method since the appearance of such discontinuities is a frequent 
and essential phenomenon of high-speed inviscid flows. 

In the following sections, the governing equations and numerical scheme are first 
summarized from [ 1 ], which is then followed by a detailed description of each of the 
test cases. 


2. Governing Equations 


The one-dimensional unsteady Euler equations of a perfect gas can be expressed 
as 


where 


+ ^ = 0, m = 1,2,3 (2.1) 

ot ox 

Ui = P, u 2 = pv, u 3 = + \pv i 2 (2.2) 

'J — 1 Z 

fi = U 2 (2.3) 

h = (7 - IK + 1(3 - 7)^ (2-4) 

r U 2 u 3 1 ,^(«2) 3 / 0 c\ 

h = "2 (7_1) W (2 ’ 5) 


with p, v,p and 7 being the density, velocity, static pressure and constant specfic heat 
ratio, respectively. 

Let x\ = x, and x 2 = t be considered as the coordinates of a two-dimensional 
Euclidean space E 2 . By using Gauss’ divergence theorem in the space-time E 2 , it can 
be shown that Eq. (2.1) is the differential form of the integral conservation law 


i h m ■ ds = 0, m = 1,2,3. 

JS{V) 


( 2 , 6 ) 


here (i) S(V ) is the boundary of an arbitrary space-time region V in E 2 , (ii) ds = dan 
with da and n, respectively, being the area and the outward unit normal of a surface 
element on 5(V), and (iii) h m = m = 1,2,3, are the space-time mass, 


3 



momentum, and energy current density vectors, respectively. Note that (i) h m • ds is 
the space-time flux of h m leaving the region V through the surface element ds, and 
(ii) all mathematical operations can be carried out as though E 2 were an ordinary 
two-dimensional Euclidean space. 


3. Euler Solver 


3(a). Marching Schemes for Interior Mesh Points 

Let H denote the set of mesh points (j, n) in E2 (dots in Fig. 1(a)) where n = 
0,±l/2, ±1, ±3/2, ±2, ±5/2, • • •, and, for each n, j = n ± 1/2, n ± 3/2, n ± 5/2, • • •. 
There is a solution element (SE) associated with each (j,n) € Cl. Let the solution 
element SE (j, n ) be the interior of the space-time region bounded by a dashed curve 
depicted in Fig. 1(b). It includes a horizontal line segment, a vertical line segment, 
and their immediate neighborhood. For the following discussions, the exact size of 
the neighborhood does not matter. 

For any ( x,t ) € SE(/, n), u m (x,t), f m (x,t), and h m (x,t ) are approximated by 
u* m {x,t\j,n), fn(x,t]j, n), and h^(x,t;j, n), respectively. They will be defined 
shortly. Let 


<,( x , t \h n ) = (“m)" +(U„„)”(Z -I,) + (U m ,)”(( -<”) 


(3.1) 


where (u m )"? and (u m t)j are constants in SE(/, n). They can be considered as 

the 
Let 


the numerical analogues of the values of and — ^ at (xj,t n ), respectively. 

ox ot 


fm,k » 75^, m, fc = 1,2, 3 

OUk 


(3.2) 


Let (/m)” and (/ m ,Jfc)j denote the values of f m and f m ,k, respectively, when u m , m = 
1,2,3, respectively, assume the values of (u m )", m = 1,2,3. Let 


and 


(A. .)? = E (UkTjMi 

k = 1 


(/ml)" = E(/-^)*(“«)i 

A:=l 


(3,3) 


(3.4) 
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It is explained in [1] that (f mx )j and {fmt)] can be considered as the numerical 


df, 


df. 


analogues of the values of ~q~ and ~£T ( x i^ n )i respectively. As a result, it is 


dt 


assumed that 


£(*,«; = (f m )] + (f mx )](x - xj) + (f mt )](t - t n ) 


(3.5) 


Because h m 


( f m , u m ) , it is also assumed that 


K(x,t;j,n) = (fjx,t-,j,n),u'jx,t;j,nj) 


(3.6) 


Note that, by their definitions, (i) (/ m )" and (f m ,k)], rn = 1,2,3, are functions of 
(u m )], m = 1,2,3, (ii) m = 1,2,3, are functions of (u m )] and (u mx )], m = 

1,2,3, and (iii) (/ m t)” are functions of (u m )] and (u mt )], rn = 1,2,3. 

Moreover, it is assumed that, for any (x,t) € SE(j,n), u m = u^(x,t;j, n) and 
fm = fm(x,t;j,n ) satisfy (2.1), i.e., 


dfa(x,t]j,n) = 
dt dx 


(3.7) 


According to (3.1) and (3.5), Eq. (3.7) is equivalent to 


(«■*)? = -(U" (3.8) 

Because (/ mx )j are functions of (u m )” and (u mi )” , Eq. (3.8) implies that (u mt )] are 
also functions of (u m )” and (u mx )j. From this result and the facts stated following 
Eq. (3.6), one concludes that (fm)], ( fmx )" and (/ OT <)” are functions of (u m )] 

and ( u mx )]. As a result, (u m )J and (u mi )?, m = 1,2,3, are the only independent 
discrete variables needed to be solved in the current marching scheme. 

The main marching scheme developed in [1] is defined by 

(u m )f = \ [ («„)“:$ + («»)£$ + (*.)£$ - (*»)”;$ ] (3.9) 

and 

(«~)J = (O* + (2« - l)(«f«~)7 (3.10) 
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where e is a parameter that controls numerical diffusion. Also, for all (j, n ) € fi: 

(3.11) 


(*-)’ = + £(/«>* + 

vn-l/2 


/ \n— i/z / in- 1/2 

(*-» = s [C«— )&g + (<-.ag] - H-),-/. (3. 


(u c ) n = 

V mi/j 


K»)"+l/2 ~ «)?— 1/2 


Ax 


with 


( U m)j±l/2 — ( U m)j± 1/2 + 2 (“ mt )j±l/2 


n— 1/2 


Let (i) 


( U m)j±l/2 ( u m)j 

Ax/2 


12 ) 


(3.13) 


(3.14) 


(3.15) 


(Umx±)” — i' 

(ii) a > 0 be an adjustable parameter, and (iii) the function Wq be defined by 

W o (0,0;a) = 0 (3.16) 


and 

W 0 (x-,x + ;a) = * j*~j/ + (|ar + | 0 + |x_| a > 0) (3.17) 

Then the main scheme can be modified by replacing (u^)" with 

(<*)? = W’o[(»-x.)J,(u.« + )J;o] (3.18) 

This modification is introduced such that numerical oscillations near a discontinuity 
can be suppressed without smearing the discontinuity. It also has the advantage that 
the modification has no discernable effect on the smooth part of the solution [1]. 


3(b). Marching Scheme for Boundary Mesh Points 

The marching schemes defined in Sec. 3(a) are applicable only for interior points 
(i) n ) € H. The boundary mesh points enter into the picture when shock-tube prob- 
lems with closed ends are considered later in the current paper. As a preliminary for 
deriving the marching scheme for boundary mesh points, we shall review the concept 
of conservation element introduced in [1]. 
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Let space-time E 2 be divided into nonoverlapping rectangular regions referred to 
as conservation elements (CEs). As depicted in Figs. l(c)-(e), conservation elements 
of three different types enter the current discussion. Conservation element CE_(j, n) 
is associated with a right-boundary mesh point (j, n) € ft. Conservation element 
CE+(j, n) is associated with a left-boundary mesh point (j, n) € ft. Finally, conserva- 
tion element CE(/, n) is associated with an interior mesh point. Note that (i) the CEs 
associated with all ( j,n ) € ft fill the entire computation domain, (ii) The boundary 
of CE_(/, n ) (CE+(j, n)) is formed by the subsets of SE(j, n ) and SE(/ — 1/2, n — 1/2) 
(SE (j -f 1/2 , n — 1/2)), and (iii) the boundary of CE(j,n) is formed by the subsets of 
SE(j», SE (j - 1/2, n - 1/2) and SE(j + 1/2, n - 1/2). 

We assumed that (i) 


f 

JS(CE(j,n)) 


• ds = 0 


(3.19) 


for any interior mesh point (j, n) 6 ft, (ii) 


f 

JS(CE-(j,n)) 




ds = 0 


(3.20) 


for any right-boundary mesh point (j, n) € ft, and (iii) 


f 

JS(CE+(j,n)) 


ds = 0 


(3.21) 


for any left-boundary mesh point (j, n) € ft. Equations. (3.19)— (3.21 ) imply that, 
for m = 1,2,3, the total flux leaving the boundary of any conservation element is 
zero. Moreover, because the flux at any interface separating two neighboring CEs is 
calculated using the information from a single SE, the local conservation conditions 
Eqs. (3.19)— (3.21) will lead to a global conservation relation, i.e., the total flux leaving 
the boundary of any space-time region that is the union of any combination of CEs 
will also vanish. 

With the aid of Eqs. (4.24) and (4.25) in [1], Eqs. (3.19)-(3.21) imply, respectively, 
(i) Eq. (3.9), which is valid for an interior mesh point (j, n) € ft, (ii) 

- + £(/„)? - + («»)? = K)":$ + (3.22) 
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(3.23) 


where (j, n ) is a right-boundary mesh point € ft, and (iii) 

j( u ’“)j n ~ + 4 Ax = ( Um )j'+l /2 — ( 5 » n )j+ 1/2 

where (j,n) is a left-boundary mesh point € ft. Equation (3.9) is a part of the 
marching schemes defined in Sec. 3(a). Note that derivation of Eq. (3.10) or its 
modified version involves considerations other than flux conservation. The marching 
schemes for boundary mesh points will be derived using Eqs. (3.22) and (3.23), and 
the appropriate boundary conditions to be imposed. 

It is assumed that, for all t > 0, 

* = 0 ' 1 = °' l = ° (324) 

at a closed end. Using Eqs. (2.1)— (2.5) and (3.24), one has 


y-2 = 0, /i = h - 0, h - (7 - 1)“3 
dfi df z 


J A v J & |"j 

~~Err ~ “oT ~ u ’ 


dt dt 


d /2 f ,,u 3 dui 
= -7(7 - 1)- 


and 


dt 

dui du 3 


dx 


= 0 


(3.25) 

(3.26) 

(3.27) 


dx dx 

at a closed end. It is assumed that the numerical analogues of Eqs. (3.25)-(3.27) are 
also valid, i.e., 


m; = o, (A)" = (Mj = 0, (/,); = (7 -i)(« 3 )j (3.28) 

(/..)” = (A)* = 0. (Mi = -7 (l - l)(>o/«i)’(«to)* (3-29) 

and 

(«i*)" = («3*)? = 0 (3.30) 

for any boundary mesh point (j, n). Here (u 3 /ui)^ = (u 3 )"/(mi)" . Combining Eqs. 
(3.22), (3.23) and (3.28)-(3.30), one obtains the marching scheme for boundary mesh 
points, i.e., (i) 



( u m) 


n— 1/2 

j- 1/2 


+ (-Sm) 


n-1/2 

i— 1/2 


0 


for m = 1,3 
for m = 2 


(3.31) 
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and 


(^mx)j — \ 


(«-)£$ + (*-)£$ - ^(7 - i)W* 

(A t) 7 , w , xn Ax 

1a 7 t(7 “ 1)( “ 3/ “ ,) ; _ X 


for m = 2 
for m = 1 , 3 


where (j, n) is a right-boundary mesh point, and (ii) 


( u m)j — \ 


( for m=l,3 


for m = 2 


(3.32) 


(3.33) 


and 


Ml = 


(«*>)";$ - Mgg + ^(7 - i)M? 

- 7 ( 7 _ 1)(U3/U X )7 + ^ 


4Ax 


0 


for m = 2 


for m = 1, 3 


where (j, n) is & left-boundary mesh point. 


(3.34) 


4. Test Problems and discussions 


To evaluate the accuracy of the numerical schemes described in the previous sec- 
tions, the numerical results of seven shock-tube problems are compared with their 
exact solutions and/or numerical solutions obtained using upwind methods. The 
shock tubes considered in these experiments are either extended to infinity or hav- 
ing one or two closed ends, / y = 1.4, and a constant value of e = 0.5 is used for all 
computations. 

Infinitely Long Shock Tubes 

Some of the shock-tube problems presented in [3] are solved again here using the 
present CE/SE Euler solver. In these problems, the shock tube is assumed to be 
long enough to avoid the reflection of discontinuities. The numerical results are com- 
pared to both the exact solutions, plotted using solid lines, and numerical solutions 
presented in [3]. 
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In the following four problems, the space region [-1,1], unless specified otherwise, 
is always subdivided into 100 intervals with Ax = 0.02, while the size of At varies 
depending on the CFL number chosen for each problem. The CFL number is defined 
as (|u| + c) max where c = \fypfp is the speed of sound. It is evaluated using the 
known exact solution. Generally the stability condition is CFL< 1 [1]. 


A. Strong Shock Waves Propagating in Opposite Directions 
The initial conditions at t=0 are: 


. %m\ I ( 0-1, °' 1 > 15 ) for x < 0 

(p,p,M) = < 

[ (0.1, 0.1, -15) for x > 0 


(4.1) 


where M = vjc. 

Two CE/SE solutions at t = lOOAt are presented in Figs. 2 and 3. Note that, 
except in the neighborhood of x = 0, these solutions agree almost perfectly with the 
exact solutions, indicating that the oscillations can be suppressed by increasing the 
value of a. 

In [3], numerical solutions to the above test problem were generated using four 
upwind schemes. The first is the AUSM (Advection Upstream Splitting Method) 
scheme. The second is the AUSMV scheme, a flux- Vector-splitting-biased modifi- 
cation of the AUSM scheme. The third is the AUSMD scheme, a flux-Difference- 
splitting-biased modification of the AUSM scheme. The last is the AUSMDV 
scheme which is a mixture of the AUSMV and AUSMD schemes. The values of 
T (temperature) , p and p generated using the above four methods ( t = lOOAt, Ax = 
0.02 and CFL= 1.0), respectively, are plotted against x in Figs. 2(a)-(d) of [3]. Note 
that T, p and p are related by the perfect gas law. 

The AUSM solution is the worst. Its values of T fluctuate violently throughout 
the interval bounded by two shock discontinuities. The maximum deviation (= 12) 
from the exact solution occurs at x = 0. There are also serious overshoots (of T 
and p ) near shock discontinuities. The AUSMD solution is only slightly better. The 
AUSMV and AUSMDV solutions are much better. Nevertheless, both have a sharp 
temperature overshoot at x = 0. From the above observations and other comparisons, 
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one may conclude that the CE/SE method is (i) more accurate in shock resolution, 
and (ii) more capable of preventing a large temperature overshoot at x = 0, than any 
one of the above four upwind schemes. 

B. Slowly Moving Contact Surface 
The initial conditions at t=0 are: 

[ (0.125, 1.0, 0.3 xvOi) for x < 0 
(/>,p,v) = < (4.2) 

( (10.0, 1.0, 0.3 X ^/0T4) for X > 0 

Note that c = \fjpfp = \/0.14 for x > 0. 

The CE/SE solution at t = 130A2 is depicted in Fig. 4. Two solutions at t = 
130 At for the same test problem obtained using the AUSMD and AUSMV schemes 
(CFL= 0.75, and Ax = 1/15) are depicted in Figs. 3(a) and 3(b) of [3], respectively. 
The AUSMD solution is comparable to the CE/SE solution in accuracy. However, the 
values of v in the AUSMV solution has a serious oscillation problem in the interval 
(- 1 , 0 ). 


C. Sod’s Problem 


Sod’s problem is solved again in order to compare the CE/SE solution with other 
solutions presented in [3]. At t = 0, we have 


{p,P,v) 


( 1 . 0 , 1 . 0 , 0 ) 

(0.125, 0.1, 0) 


for x < 0 
for x > 0 


(4.3) 


The CE/SE solution at t = 50 At is plotted in Fig. 5. The solution at t = 50 At 
for the same test problem obtained using the AUSMDV, AUSM and Roe schemes 
(CFL=1.0 and Ax = 0.02) are presented in Figs. ll(a)-(c) of [3], respectively. For 
each of the above three schemes, it requires about 15 mesh intervals to resolve the 
contact discontinuity and 2-3 mesh intervals to resolve the shock discontinuity. On 
the other hand, the CE/SE method requires only 4-5 mesh intervals and 1 mesh 
interval, respectively, to resolve the contact and shock discontinuities. Moreover, the 
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numerical errors of the CE/SE solution in the smooth-solution regions are also less 
than those associated with other three solutions. 

D. Strong Expansion Fans 

The initial conditions at t=0 are: 

( (1.0, 2.0, — 2.5x\/2^8) for x < 0 

ObP.w) = \ , (4-4) 

( (1.0, 0.5, 2.5xv^8) for x > 0 

Note that c = \jypTp = for x < 0. 

Assuming CFL=1.0 and Ax = 0.02, the numerical values of p,p and T at t = 
30At = 0.1025 for the above test problem obtained using the AUSMDV, Godunov 
and HILE schemes, respectively, are presented in Figs. 12(a)-(c) in [3]. Because the 
CE/SE solver with e — 0.5, a = 1 and Ax = 0.02 is not stable at CFL=1.0, a direct 
comparison of the CE/SE method with the above three schemes can not be made. 
Instead, assuming e = 0.5, a = 1 and Ax = 0.02, we reach for the largest At (i.e., the 
largest CFL because Ax = 0.02 and (|u| + c) max = 5.8566 are known) such that (i) the 
CE/SE solver is stable, (ii) nAt = 0.1 where n is an interger, and (iii) the numerical 
solution at t = 0.1 is such that p > 0 and T > 0 at all mesh points. Note that 
p < 10 -3 in the neighborhood of x = 0.16. As a result, a numerical undershoot could 
cause p to become negative. It is found that At = 1 /550 (corresponding to n = 55 and 
CFL=0.532) satisfies the above requirements. The corresponding solution is plotted 
in Fig. 6. 

By comparing the CE/SE solution with the Godunov solution which is depicted 
in Fig. 7, it becomes obvious that the former are far more accurate than the latter. 
The Godunov’s values of T deviate substantially from the exact values in the interval 
(—0.5, 0.5). The Godunov’s values of p also suffer the same order of error percentage- 
wise although it is difficult to see because the value of p are very small in the same 
interval. The numerical errors associated with the AUSMDV solution are similar to 
those of Godunov’s solution while the HILE scheme produces even worse results. 

In general, the comparison with both the exact and other upwind solutions in 
this part shows that the current CE/SE Euler scheme can generate highly accurate 
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shock-tube solutions. Having demonstrated the Euler solver’s robust capability of 
capturing shock and contact discontinuities as well as simulating expansion waves, 
we are proceeding in the following to extend its application to study the end effects 
on shock-tube flows. 


Test Problems Taken from Ref. 4 

In this part, more complex shock-tube problems involving the reflection and in- 
tersection of shock waves, expansion waves and contact discontinuities are solved by 
use of the present Euler solver. The following three test cases are taken from [4] 
where the exact solutions are also given. The results presented here are obtained 
with CFL=0.82 in all three cases, whereas the value of Ax is 0.1 in cases E and F 
but is changed to 0.01 in case G for a more detailed description of the flow. 


E. Shock Tube With a Closed End 


Consider a strong shock wave traveling at M=4.6 in a constant-area duct filled 
with air [4, pp. 143-144]. When the shock reaches the closed end on the right, it will 
reflect as described in Fig. 8. 

Let the closed end be located at x = 3.5, and the initial conditions at t=0 be 
given as 


(p> Pi v) = < 


(4.5) 


(24.52, 4.852, 4.321) for x < 0 

(1.0, 1.0, 0) for x > 0 

Figures 9 and 10 show respectively the results at t=0.5 when the incident shock 
moves to the right toward the closed end, and at t=1.5 after the shock wave has been 
reflected and reversed its direction . Other than a slightly lower density observed near 
the closed end in the reflected wave, the CE/SE solution has excellent agreement with 
the exact solution at both t = 0.5 and t = 1.5. 


F. Merging of Two Shock Waves 

With reference to Fig. 11, consider two shock waves propagating to the right in 
an infinitely long tube with a strong shock of M = 3.0 behind a weaker shock of 
M = 1.5 [4, pp. 131-134]. 
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Let the initial conditions at t=0 be 


(P, P, v) 


(25.41, 7.189, 3.842) 
< (2.460, 1.863, 0.821) 
( 1 . 0 , 1 . 0 , 0 . 0 ) 


for x < —1.5 
for — 1.5 < x < 2.6 
for x > 2.6 


(4.6) 


The computed results at t=1.0, 1.25 and 1.75 are shown in Figs. 12-14 in favorable 
comparison with the exact solution. These two shock waves remain separated before 
t=l .25. At t=1.25, the faster shock catches up with the slower shock ahead and 
merging of the two waves begins. As a result of merging, a shock wave stronger than 
the original faster shock is created and is propagating to the right at a higher speed, 
leaving behind it a contact surface and a left-moving expansion fan wave, as shown 
in the plot of Fig. 15 at t=1.75. 


G. Shock Tube With Two Closed Ends 

This shock tube problem depicted in Fig. 15 constitutes of reflection of shock wave 
and expansion wave and the intersection of the reflected waves and contact surface. 
The tube has closed ends at x=0 and x=l, respectively. At t=0, a diaphragm located 
at x=0.25 separating two gases at different conditions is bursted [4, pp. 205-208]. 

With the initial conditions that 

f (20.0, 20.0, 0.0) for x < 0.25 
(P, />,*>) = < (4-7) 

[ (1.0, 1.0, 0.0) for x > 0.25 

the break of the diaphragm creates a shock wave and an expansion wave separated 
by a contact surface . 

The calculation is carried out for a period that is long enough to show successive 
reflection and intersection of shock wave, expansion wave and contact surface. When 
t= 0.09, waves have not yet reached the ends as shown in Fig. 15. When t=0.3, 
the left-moving expansion waves have been reflected from the left end but the right- 
moving shock wave has not . At t=0.4, the shock wave has already completed its 
reflection from the right end and the reflected shock wave is moving to the left. This 
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reflected shock is intersecting with the right-moving contact surface to create a left- 
moving transmitted shock, a right-moving contact surface and a right-moving shock 
wave which, when reaching the right end, is reflected again. At t=0.585, there exist in 
the flow field two shock waves, one contact surface and a region of expansion waves. 
Numerical results plotted in Figs. 16—19 correctly display the aforementioned flow 
phenomena at those four time instants. 

The excellent agreement in all cases demonstrates the ability of the CE/SE scheme 
to accurately capture weak and strong waves and discontinuities even in the presence 
of their reflections and intersections. It again verifies the robustness of this scheme. 

5. Conclusion 

The Space-Time Conservation Element and Solution Element method has been 
applied to solve the Euler equations. Seven numerical examples on one-dimensional 
shock-tube flows have been used to demonstrate not only the advantages of simplic- 
ity, accuracy and generality of this method, but also its ability of achieving high 
resolution of discontinuities even in the presence of wave intersections and reflections 
without any ad hoc techniques. It shows that the CE/SE method is a powerful tool 
of high potential for solving general fluid dynamic problems. Under investigation is 
the extension of this method to tackle flow problems of higher dimensions. 
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Fig. 1 The SEs and CEs of the CE/SE method 
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Fig. 2 The CE/SE solution for case A at t=100At 

(CFL=1.0, AX =0.02, a=l) 
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Fig. 3 The CE/SE solution for case A at t=100At 

(CFL=1.0, AX =0.02, a =3) 
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Fig, 4 The CE/SE solution for case B at t=130At 

(CFL=0.75, AX=1/15, a=l) 
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Fig. 5 The CE/SE solution for case C at t=50At 
(CFL= 1.0, AX =0.02, a=l) 
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Fig. 6 The CE/SE solution for case D at t=0.1 
(CFL= 0.532, AX =0.02, o=l) 
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Fig. 7 The Godunov solution for case D at t=0.1025 

(CFL=1.0, AX=0.02) 
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Fig. 8 Shock reflection from a closed end 
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Pressure Velocity Density 


Fig. 9 The CE/SE solution for case E at t=0.5 
(CFL=0.82, AX=0.1, a=3) 
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Fig. 10 The CE/SE solution for case E at t=1.5 
(CFL=0.B2, AX=0.1, a =3) 
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Fig. 11 Merging of two shock waves 
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Fig. 12 The CE/SE solution for case F at t=1.0 

(CFL=0.82, AX=0.1, o=3) 




28 













Fig. 13 The CE/SE solution for case F at t=1.25 

(CFL=0.82, AX =0.1, a=3) 
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Fig. 14 The CE/SE solution for case F at t=1.75 

(CFL=0.82, AX=0.1, a=3) 
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Fig. 15 Shock tube with two closed ends 
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